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ABSTRAGT 

The Apache Point Observatory Galactic Evolution Experiment (APOGEE), part of the Sloan Digital 
Sky Survey HI, explores the stellar populations of the Milky Way using the Sloan 2.5-m telescope 
linked to a high resolution (R^22,500), near-infrared (1.51-1.70 fim) spectrograph with 300 optical 
fibers. Eor over 150,000 predominantly red giant branch stars that APOGEE targeted across the 
Galactic bulge, disks and halo, the collected high S/N (>100 per half-resolution element) spectra 
provide accurate (^0.1 km s“^) radial velocities, stellar atmospheric parameters, and precise (<0.1 
dex) chemical abundances for about 15 chemical species. Here we describe the basic APOGEE data 
reduction software that reduces multiple 3D raw data cubes into calibrated, well-sampled, combined 
ID spectra, as implemented for the SDSS-HI/APOGEE data releases (DRIO, DRll and DR12). The 
processing of the near-IR spectral data of APOGEE presents some challenges for reduction, including 
automated sky subtraction and telluric correction over a 3°-diameter field and the combination of 
spectrally dithered spectra. We also discuss areas for future improvement. 

Subject headings: methods: data analysis - techniques: image processing - Galaxies: kinematics and 
dynamics - Galaxies: Local Group - Galaxy abundances - Galaxy: halo - stars: 
abundances 


1. INTRODUCTION 

The SDSS-HI/APOGEE project obtained high resolu¬ 
tion IR spectra of over 150,000 Milky Way stars during 
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the period 2011-2014 , as des cribed in Alam et al. (2015) 
and Majewski et al. (2015). It is based on a bench- 
mounted spectrograph operating in cryogenic conditions 
that can obtain 300 simultaneous spectra covering wave¬ 
lengths of 1.51-1 .70/im for light fed to it from the Sloan 
2.5-m telescope (|Gunn et al. 2006) via 40-m long opti¬ 
cal fibers. The spectra are recorded onto three separate 
Hawaii-2RG (H2RG) arrays, where each array covers a 
different wavelength span, with small gaps in wavelength 
coverage between. The instrum ental design and per for- 
mance are described in detail in Wilson et al. (2015). 

Given the large amount of data expected for the entire 
survey, development of an automated reduction pipeline 
was essential. However, reduction of APOGEE spectra 
is not straightforward for several reasons: 


• In the near-IR, telluric absorption is significant and 
both spatially and temporally variable; telluric fea¬ 
tures affect a significant fraction of the APOGEE 
observed spectrum. 


• The sky brightness, dominated by OH lines, is tem¬ 
porally and spatially variable. 

• The instrument design delivers slightly undersam¬ 
pled spectra at its short wavelength end. To avoid 
issues with undersampling for the stellar parame¬ 
ter and abundance determinations, data are taken 
at two different dither positions, where the entire 
detector assembly is shifted by ^0.5 pixel between. 

• The Teledyne H2RG arrays that are in the 
APOGEE instrument have some performance com¬ 
plications; in particular, some regions show signif¬ 
icant persistence, where previous exposure to light 
affects the subsequent behavior of the detector. 
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This paper describes the state of the APOGEE data 
reduction pipeline as it has been used to produce data 
contained in the data releases for SDSS-III APOGEE, 
DRIQ (June 2013: Ahn et aL]|2Q14 ) and DR12 (January 


2015; Alam et ah 2015| ). As described below, there are 


still some areas for potential improvement in the pipeline, 
but the main goal of this paper is to document the meth¬ 
ods used to process the APOGEE database for its public 


are discussed in Holtzman et al. 

( 

2015| 

tails on the stellar parameter and 

.abund 

presented in Garcfa-Perez et al. 

(2015) 


cussed in Holtzman et ah (2015) for DR12 and |Me 


et al.lpOl^ for DRlOp] 


leszaros 


I'he layout of this paper is as follows. Section [^details 
the survey operations and data taking. An overview of 
the reduction pipeline is given in Section Sections [4}|^ 
describe the three steps of APRED that reduce observa¬ 
tions of an individual plate on a single night: (1) AP3D 
reduces the 3D data cubes to 2D images, (2) AP2D ex¬ 
tracts the 300 spectra and determines wavelength zero- 
points, and (3) APIDVISIT performs sky corrections, 
combines multiple dither exposures, and determines ini¬ 
tial radial velocities. APSTAR combines spectra of indi¬ 
vidual stars on the rest-frame and determines accurate 
radial velocities (Section!^. Radial velocity determina¬ 
tion for both the APlDvTSIT and APSTAR steps are 
described in detail in Sectio n Isl Access to the data prod¬ 
ucts is described in Section |^and, finally, a summary is 
given in Section Eigure]^ shows a flowchart of the 
pipeline processing steps. 

2. SURVEY OPERATIONS AND DATA TAKING 

Initial commissioning data were taken starting in April 
2011. The initial data showed some issues with the inter¬ 
nal optics and with relative focus of the three detectors, 
so the instrument was opened during summer 2011 to 
improve these issues. It was subsequently cooled in Au¬ 
gust 2011, and official survey-quality data began to be 
collected after this time. The instrument remained cold, 
and in the same optical configuration over the course of 
the entire survey and is therefore quite uniform. 

To collect data, the 300 APOGEE optical fibers are 
coupled to standard Sloan 2.5-m plug plates. Eor rou¬ 
tine first year science observations, 230 of the fibers are 
placed on science targets (almost all stars); 35 additional 
fibers are placed on blue stars to be used to measure tel¬ 
luric absorption, and 35 fibers are placed on sky regions 
without objects. Targeting is based almost e ntirely on 
the 2MASS catalog and is described in detail in |ZasowsE 
et al. (2013). Eor most survey fields, the fibers are dis- 
tributed over a 3° (diameter) held of view, although for 
some low declination, high airmass fields, a smaller held 
is used to minimize differential refraction effects. 

Data are collected from the SDSS telescope using 
the standard SDSS telescope/instrument interface, STUI 
(SDSS Telescope User Interface). The normal mode of 
operation is to take individual exposures of 500s dura¬ 
tion. Exposures are taken at either of two dither posi¬ 
tions (“A” or “B”), where the detectors are nominally 

DRII was an internal collaboration data release but followed 
the same procedures. 


moved by ^0.5 pixels between exposures (although early 
survey data had a somewhat smaller shift of ^0.4 due 
to technical issues but this did not adversely affect the 
reduction). The standard observing sequence collects 
2 ABBA exposure sequences per plate, which leads to 
slightly over one hour of exposur e per plate on a night . 
Eor most fields (but not all, see Zasowski et al. 2013), 
the target exposure time is three hours, but these are 
collected over three different visits to the held, spread 
out in time, to enable identification of radial velocity 
variation arising from stellar binarity. Individual visits 
to a held are identified by a plate identification and an 
MJEp^ Individual visits to an object are identified by a 
plate, MJD, and fiber number; a given star will not gen¬ 
erally be observed in the same fiber in subsequent visits, 
since plates are typically replugged between visits. 

The data collection system continually reads the de¬ 
tectors non-destructively as the charge is being accumu¬ 
lated, at slightly more than 10 seconds between reads, so 
the 500s exposures are composed of a series of 47 read¬ 
outs. Since these “up-the-ramp” readouts are accessi¬ 
ble as the exposure is proceeding, it is possible to ana¬ 
lyze count rates as the exposure is accumulating. This 
is done by “quicklook” software that communicates this 
information to the SDSS observers. To date, this has 
largely been used for informational monitoring, and total 
exposure times have been fixed to 500s. Under poor con¬ 
ditions, a third ABBA sequence is sometimes obtained. 
During the final visit (usually the third) to a held, only 
one ABBA sequence is taken if it is likely to be sufficient 
to reach the total desired accumulated S/N. 

After exposures are finished, a quick reduction is done 
to provide observers with some roughly reduced data 
(i.e., extracted ID spectra) to inspect. The quick re¬ 
duction software also takes the files with the individual 
readouts, bundles them into three data cubes (one for 
each detector), and compresses them (see next section). 
Additionally, information about the exposure and quick- 
reduced spectra are inserted into a database running on 
the mountain. This database is used to monitor progress 
of the observations on each held, and is the basis for the 
autoscheduler, which determines the plan for plugging 
and observing of new fields. A web application allows 
for a graphical interface to this database. 

Galibration data are obtained by coupling the fiber 
bundle from the instrument to a fiber bundle that leads 
to an integrating sphere with calibration sources. Three 
calibration sources are available: a continuum source, 
and Thorium-Argon-Neon (ThArNe) and Uraninum- 
Neon (UNe) hollow-cathode lamps. In addition, several 
IR LEDs were installed on a cold shutter mechanism that 
was installed in summer 2011. These are located down¬ 
stream of the internal slithead, and provide roughly uni¬ 
form illumination of the detectors that can be used to 
determine pixel-to-pixel sensitivity variations; we refer 
to such frames as internal flats. 

Some calibration data are taken on a daily basis, 
mostly for instrument performance monitoring. At the 
end of the afternoon, a few test frames of a continuum 
source and some line lamp exposures are taken for quick 


Modified Julian Date (MJD) = Julian Date (JD) — 2400000.5. 
The MJD used by SDSS-III is MJD+0.3 days so that the “day” 
increments in the afternoon at Apache Point Observatory. 
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Fig. 1.— Flowchart for the APOGEE data reduction pipeline listing the steps of the main stages and the resulting data products. 


inspection by the observers to confirm routine instrument 
performance. At the end of each night, a more complete 
set of calibration data are obtained; this set includes sev¬ 
eral long dark frames, lamp exposures (ThArNe and UNe 
at both dither positions), and several internal fiat fields. 
More extended sequences of calibration data were taken 
near the beginning of the survey and were repeated pe¬ 
riodically throughout the duration of the survey. 

In addition, on each observing night, 4 exposures (1 
X ABBA) are taken of a random pointing on the sky 
so that all fibers are illuminated by sky (predominantly 
OH lines). These frames are used to characterize and 
monitor the image quality and the related spectral line 
spread function (LSE). 

2.1. Data volume and compression 

The nightly volume of data collected is significant. 
Each standard exposure has 47 readouts of three 
2048x2048 chips, and the instrument computer also col¬ 
lects an additional 2048x2048 array that provides bias 
information. This leads to roughly 1.5 GB per exposure. 
Eor multiple exposures and multiple plates, plus associ¬ 
ated calibration data, this leads to of order ^100 GB of 
data per full night of APOGEE observing. 

All of the raw up-the-ramp APOGEE data are kept 
and transferred daily off the mountain to the SAS. To 
speed up the data transfer off the mountain and reduce 
overall disk space the raw data are compressed using a 


custom designed algorithm. This algorithm takes advan¬ 
tage of the fact that successive reads of the arrays are 
very similar; as a result, the sequence of difference im¬ 
ages has relatively smaller dynamic range and can be 
compressed efficiently. Three steps are used to compress 
the up-the-ramp data cubes: 


1. The detector reads are converted into difference im¬ 
ages resulting in N^eads integer images (Nreads-1 dif¬ 
ference images and the first read). 

2. The average difference image (rounded to integers) 
is computed and subtracted from the difference im¬ 
ages resulting in Nreads+1 integer images (Nreads- 
1 “residual” images, one average difference image, 
and the first read). The Nreads+1 integer images 
are written to a multi-exension EITS file. 


3. The EIT S file is compresse d using the EPAGKp^ 


routines (Pence et al. 
compression algorithniT 


2010) and the lossless Rice 


These custom APOGEE compressed files (a separate one 
for each of the three arrays) are saved to disk with “.apz” 
extensions and are on average compressed by a factor of 
^2. The theoretical best compression rate of data of our 
noise-level (^20) and bits per pixel (16) is ^2.6 (Pence 


I http://heasarc.nasa.gov/fitsio/fpack/ 
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et al. 2009), although in practice the best algorithms 
will obtains compressions of ^2.2. Therefore, our cus¬ 
tom compression algorithm is about as good as can be 
expected for our data. 

Compression and uncompression algorithms are imple¬ 
mented in the APZIP and APUNZIP custom IDL pro¬ 
cedures. 


3. PIPELINE OVERVIEW 

Data reduction is run off-site on the raw, compressed 
data downloaded from the SAS. There are two main 
stages for basic data reduction: 

1. APRED reduces observations of an individual plate 
on an individual night in three steps: 

(a) AP3D reduces the data cubes to 2D images, 
applying detector calibration products in the 
process. A separate error image and a mask 
of bad pixels are also calculated. 

(b) AP2D extracts the spectra from each 2D im¬ 
age to produce 300 well-sampled ID spectra 
and corrects them for throughput variations. 
AP2D also performs wavelength calibration 
using exposure-specific wavelength zeropoints 
determined from the positions of sky lines and 
wavelength solutions previously derived from 
emission-line lamp exposures. Error images 
and masks are also produced. 

(c) APIDVISIT measures accurate dither shifts 
between exposures in a visit, corrects individ¬ 
ual exposures for sky emission and absorption, 
and combines multiple dithered exposures in 
a visit to produce 300 ID spectra. The code 
then determines an initial radial velocity es¬ 
timate for each stellar object using a best¬ 
matching stellar template. Outputs include 
sky-corrected spectra, as well as pixel-by-pixel 
errors, mask information, and wavelength ar¬ 
ray, and the correction spectra used for the 
sky correction. In addition, information about 
the observed dithering and dither combination 
are output. 


2. APSTAR is run after multiple individual visits 
have been obtained. It resamples the spectra onto 
a fixed wavelength grid (constant dispersion in 
log(A)), correcting each visit for the visit-specific 
radial velocity, and coadds the spectra. Using 
the combined spectra, the code derives relative ra¬ 
dial velocities via cross-correlation of each individ¬ 
ual spectrum with the combined spectrum, and 
puts them on an absolute scale by cross-correlating 
the combined spectrum against a best-matching 
template spectrum. As a check, individual visit 
RVs are rederived by cross-correlating each visit 
spectrum against the common template that best 
matches the combined spectrum. 


A flowchart of these processing steps and the resulting 
data products are shown in Eigure[2 
A third pipeline stage determines stellar parameters 
and chemical abundances in the APOGEE Stellar Pa¬ 
rameters and C hemical Abundances pipeli ne (ASPCAP), 
as described in Garcfa-Perez et al. (2015). 


The software used for the APOGEE data reduction 
pipeline is almost exclusively implemented in the In¬ 
teractive Data Language (IDL) The code is archived 
and managed through use of the SDSS-III software 
repository using the software management package Sub¬ 
version (SVN). While it was not designed for gen¬ 
eral public usage, the DR12 version of the code is 
publicly available http://www.sdss3.org/svn/repo/ 
apogee/apogeereduce/, 

4. AP3D: REDUCTION OF DATA CUBE TO 2D IMAGE 

During an exposure the three APOGEE arrays are read 
out non-destructively every ^10.6 seconds in sample-up- 
the-ramp (SUTR) mode. In addition to providing the 
opportunity for inspection of data as it is being accumu¬ 
lated, the SUTR can be used advantageously to: 

• reduce the read noise by using multiple measure- 
ments of the electron s as they are accumulated 
( ^anscher et al.||2QQ7 ); 

• detect and correct cosmic rays; 

• potentially correct saturated pixels if there are 
enough (^3) unsaturated reads to measure the flux 
rate, and the flux rate is assumed constant in time. 

As discussed above, the separate APOGEE raw read¬ 
outs are stored in a data cube (one per array) and subse¬ 
quently compressed. In AP3D, these data cubes are “col¬ 
lapsed” into 2D images. Basic calibration is also done at 
this stage, leading to these main steps: 

1. reference pixel correction, 

2. linearity correction (currently not implemented), 

3. dark subtraction, 

4. cosmic ray detection and repair; saturated pixel 
correction (currently not implemented) 

5. collapse to 2D image, 

6. flat helding, 

7. construction of error array and bad pixel mask. 

These are ah described in more detail below. 

4.1. Detector electronics and reference pixel correction 

Each of the three Teledyne Hawaii 2RG detectors 
(each with 2048 x 2048 18/im pixels) are read in parallel 
through 4 different channels per chip, with each “quad¬ 
rant” being 512 x 2048 in size. 

The voltage bias for the Hawaii 2RG arrays can drift 
slowly over time, but reference pixels have been imple¬ 
ment to correct for this effect. There are two types of 
reference pixels: (1) a perimeter of 4 pixels around each 
array (“embedded” reference pixels) that are not “ac¬ 
tive” but are read out the same way as the rest of the 
array (via 4 output channels per array). (2) A single 
reference pixel for each array that is read out with its 
own readout port, and is called the “reference output”. 

A product of Exelis Visual Information Solutions formerly 
ITT Visual Information Systems and Research Systems, Inc. 
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This output channel is read out in parallel and at the 
same rate as the other four output channels (leading to 
five altogether) producing a separate 512 x 2048 image. 
In the raw data cubes this extra 512 x 2048 “reference” 
image is attached to the end of the regular read image 
to produce a 2560 x 2048 array. The reference image is 
useful to correct for electronic “ghosts” that are created 
when very high counts from a single output affect the 
other three. 

Both types of reference pixels are used in the APOGEE 
reduction. Eirst, the 512 x 2048 reference arrays are sub¬ 
tracted from each quadrant. Second, vertical ramps are 
created for each quadrant using means of the bottom/top 
reference pixels and then subtracted. Einally, horizon¬ 
tal ramps are created for the entire array using 50-pixel 
smoothed values from the left/right reference pixels and 
subtracted. This is performed separately for each read¬ 
out and detector. 


4.2. Linearity 

Most infrared detectors have some degree of nonlinear¬ 
ity, with th e sensitivity changi ng slightly as charge is ac¬ 
cumulated (Kubik et al.||2014 ). Linearity corrections are 
likely to be less important for our spectroscopic analysis 
because the dynamic range of a given spectrum, espe¬ 
cially over the portion of any individual spectral features, 
is generally relatively small, and we have no requirement 
for high accuracy in relative flux between differ ent ob¬ 
jects. Takin g the nonlinearity coefficients from |Kubik| 
et al.| ( |2Q14 ), the error in the relative depth of a 20% 


(of the continuum) absorption line for a spectrum with 
a continuum of 10,000 ADU is ^0.4%. The large ma¬ 
jority of our stellar spectra have continua (in individual 
exposures) below this level, and, therefore, the effects of 
nonlinearity are minor. 

We have made some initial tests for nonlinearity using 
internal flat field data cubes, under the assumption that 
the LED light source is stable (which it appears to be, 
judging from the repeatability of light levels in succes¬ 
sive exposures). These data suggest that there may be 
some small level of nonlinearity, but characterizing it is 
complicated by the behavior of some pixels at low light 
levels. As a result, we have chosen not to implement any 
nonlinearity correction at this time. 

Additionally, as discussed below, some regions of two of 
the detectors suffer from a significant persistence effect, 
where the amount of charge deposited can be affected by 
the previous exposure. This effect is significantly larger 
than any expected non-linearities in these regions. 

Gonsequently, although the pipeline has an implemen¬ 
tation for a linearity correction, we have not applied such 
a correction for the DR10-DR12 data. 


4.3. Dark current 

The dark current is derived from multiple 60 read ex¬ 
posures with the internal cold shutter closed. Eigure 
shows the distribution of the dark current rate (counts 
per read) of all of the pixels. Three panels show the his¬ 
togram of dark rates for each of the three chips, whereas 
the lower right shows the cumulative distribution across 
all pixels for all three chips. While most pixels have a 
dark rate below 0.5 counts/read, there is a tail up to 
high dark rates for some pixels. In fact, the typical dark 



0.01 0.10 1.00 10.00100.00 
dork per read 



0.01 0.10 1.00 10.00100.00 
dork per reod 



0.01 0.10 1.00 10.00100.00 
dark per read 



0.01 0.10 1.00 10.00100.00 
dark per read 


Fig. 2. — Histogram of dark current rates (DN/read). Three 
panels show histograms for the three chips respectively; the lower 
right shows the cumulative histograms for all three. Data are shown 
for two different dates separated by over two years. 


current is significantly lower than 0.5 counts per read be¬ 
cause, even after averaging 20-30 frames to reduce the 
noise, readout noise dominates over dark current at this 
level. The middle array has a section of higher dark 
current than the other two, and this is reflected in the 



adjacent pixels. 

To correct for dark current, “superdark” frames are 
constructed by taking the median of 20 long dark frames. 
To allow for the possibility that dark current may not ac¬ 
cumulate linearly with time, the superdark is constructed 
for each up-the-ramp readout, so the superdark calibra¬ 
tion frame is a data cube, with the appropriate slice sub¬ 
tracted from the corresponding readout of each science 
frames. 

Eor the current analysis, we have used a single super¬ 
dark, constructed from data taken on MJD 56118. Eigure 

t (top and lower-left panels) shows the histograms from 
is date in black, as well as the histograms from darks 
taken at the end of the SDSS-III/APOGEE survey, on 
MJD 56853 (in red). Gomparison of these histograms, as 
well as direct comparison of the dark images, show that 
the dark current is quite stable, with only a few pixels 
changing their dark rate significantly. Daily dark frames 
are taken as part of the normal calibration, and future 
analysis may use these to implement correction for small 
changes in dark rates and the number of hot pixels. 

All pixels with a dark rate exceeding 10 count s/read 
are marked as bad pixels. Neighbors of any of these are 
marked as bad if their dark rates exceed 2.5 counts/read. 

4.4. Cosmic Ray and Saturated Pixel Correction 
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Row 488 Column 87 





Fig. 3. — Example of cosmic ray detection and “repair”. (Top) 
Counts and (bottom) count rate as a function of read number. 
The original values (red), repaired values (blue), detected cosmic 
ray (green), and linear fit to the fixed values (black) are shown. 


One advantage of using SUTR sampling is that cosmic 
rays can be detected and removed. Each pixel is searched 
for positive jumps that could correspond to cosmic rays. 
Eirst, the array of read values is turned into difference 
counts (between successive reads). Second, a median fil¬ 
tered version of the difference counts is created (using 
a median filter of 11 reads) that can be used to remove 
any flux rate variations over time (e.g., from seeing varia¬ 
tions or clouds). Then, a “local” scatter of the difference 
counts around the local median value is measured using 
a robust standard deviation. Any difference count values 
larger than 10 x the local scatter above the local median 
(and 10 X above the noise level) are flagged as cosmic 
rays. The cosmic rays are corrected by replacing their 
difference values with the local median for that pixel. 
Einally, the counts array is reconstructed by adding the 
cumulatively summed difference counts to the first read 
value. The pixels with detected cosmic rays are flagged in 
the mask image. This cosmic ray detection method will 
miss some very weak cosmic rays but should catch most 
of the large ones that are more likely to affect the data. 
Eigurej^ shows an example of the cosmic ray correction 
process. Eigureffl shows the distribution of detected cos¬ 
mic ray rates in me DR12 data for the three arrays. The 
median detected cosmic rays in object exposures (500s) 
are 430/470/470 (blue/green/red). 

Another advantage of using SUTR sampling is that 
measurable signal is recorded even for pixels that end up 
being saturated after the full exposure time. If as few 
as ^3-4 reads are unsaturated, then the flux rate can be 
measured and used to extrapolate the counts to the end 
of the exposure. However, this extrapolation assumes 
that the count rate is stable, which is not the case in sub- 
optimal observing conditions. While it might be possible 
to characterize the count rate variations using other pix¬ 
els, the situation is complicated because different pixels 
may have different variations in count rate, e.g., the rate 


Fig. 4. — Histogram of cosmic rays detected per exposure time 
on the three APOGEE detectors ('^IShO mm^). The median val¬ 
ues are 0.92/1.00/0.99 (blue/green/red). The images with shorter 
exposure times have slightly higher rates likely because of lower 
Poisson noise that allows for the detection of weaker cosmic rays. 

in spectral regions including sky emission are likely to 
vary in a different way from the rate in regions of stellar 
signal. 

The pipeline currently corrects any saturated pixels as¬ 
suming a constant flux rate, but flags such pixels as hav¬ 
ing been corrected. Because of the possibility that this 
correction is a poor approximation, subsequent stages of 
the reduction treat these pixels as bad. 


4.5. Collapse to 2D Image 

As is standard with non-destructive readout IR arrays, 
the array reset at the beginning of each exposure intro¬ 
duces a pedestal value with significant “noise” (i.e., pixel- 
to-pixel variations), so multiple non-destructive readouts 
are required to remove the baseline values and this “re¬ 
set noise”. To minimize this non-negligible readout noise, 
2D IR images are typically calculated from the 3D data 
cubes using either Eowler sampling, where some number 
of readouts at the beginning and end of the exposure 
are averaged and a difference image created, or using the 
full up-the-ramp sampling to determine a mean count 
rate that, when multiplied by the exposure time, gives 
the total counts per pixel. The simplest Eowler sampling, 
using just one readout at the beginning and one at the 
end, corresponds to correlated double sampling (CDS). 
A detailed analysis of t he noise using differen t readout 
methods is presented in Rauscher et al. (2007). 

Both Eowler sampling and up-the-ramp analysis are 
implemented in the pipeline. Up-the-ramp is used for all 
of the data except for the “dome” flats; these use simple 
CDS sampling because the lamp is only turned on for a 
few seconds to accumulate the desired number of counts; 
because the count rate is thus highly non-uniform, up- 
the-ramp sampling fails in this case. In addition, the 
relatively large number of counts in these exposures is 
not significantly affected by readout noise. 
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Fig. 5.— The APOGEE DR12 flat fleld images. Note the different colorbar ranges for each chip. The standard deviations are 3.0, 8.2 
and 7.6% for the red, green and blue detectors. 
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Fig. 6. — Histogram of the ratio of two superffats, taken over two 
years apart. The red lines show gaussian fits to the histograms, 
and the a of the fit is shown. Photon statistics alone are expected 
to yield a = 0.002 — 0.003. 


4.6. Flat Fielding 

Elat fields to correct for pixel to pixel sensitivity vari¬ 
ations are constructed using a median of ten internal fiat 
fields and stored in the apElat calibration files. Eigure 
shows the observed sensitivity variations of the three 
Hawaii 2RG detectors which can be quite large. The 
middle (“green”) chip shows a rim of lower QE (quantum 
efficiency) giving the appearance of a “thumb-print”. 

The 2D images are fiat fielded using a “superfiat” con¬ 
structed by the average of 20 internal flat field frames, 
with large scale structure removed; these are stored in 
apElat files. The count levels in the individual flats are 
such that the combined flat should have S/N> 500, or 


an rms of <0.2%. 

The flat field is very stable in time. This is demon¬ 
strated in Eigure which shows histograms of the ra¬ 
tio of a superfiat constructed on MJD 56037 with those 
taken at the end of the survey on MJD 56852. Given 
the S/N of the individual superfiats, the ratio of the two 
flats is expected to have a ^ 0.003 on the basis of pho¬ 
ton statistics only. Eigure shows that the distribution 
of values in the ratio closely matches a gaussian with 
width comparable to this, demonstrating that the flats 
are likely to be extremely stable. We note some devia¬ 
tions in the ’b’ chip; upon inspection, these arise from 
the rim noted above, which is likely more re nte d to per¬ 
sistence effects in this detector (see Section 4.8) than to 
temporal QE variations. 

Any pixels in the superflat that have sensitivity less 
than 75% of pixels in their vicinity are marked as bad 
pixels, and saved in the bad pixel masks. 


4.7. Gain, Noise Model and Bad Pixel Mask 

The detectors are read in parallel through 4 different 
channels per chip. The inverse gain (e“ / DN) for each 
channel in each chip is derived from pairs of internal flat 
fields using the measured variance. In 50 different in¬ 
tensity bins with different mean intensities (/), the vari¬ 
ance in the difference image is measured and the inverse 
gain calculated from 2//cr^. We use a gain of 1.9 for all 
quadrants. This gain calculation will be updated in the 
future. 

A noise model is constructed by combining in quadra¬ 
ture the Poisson noise of the 2D image, the Poisson noise 
of the dark ima ge, and the sampling r ead noise (using the 
equations from Rauscher et al. 2007). This noise image 
is placed in the “error array”. 

Bad pixel mask files for each detector are created using 
pixels marked as bad during the construction of the su¬ 
perdark and superflat frames. The bad pixel masks are 
saved as bitmasks to preserve the reason that any given 
pixel was marked as bad, and these values are propagated 
via the bitmasks to the reduced data frames. 

The bad pixel mask from the apBPM calibration prod¬ 
uct is used to mask out (set to NAN) bad pixels from each 
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Median Continuum of Sky Fibers in Blue Detector 




H-magnitude Difference (AH) 


Fig. 7.— The median continuum of sky fibers on the blue detec¬ 
tor versus exposure number in a visit (median across 1638 survey 
visits) for moderately high persistence region (green), highest per¬ 
sistence region (red), and normal persistence region (blue). The 
magnitude and a rate of (temporal) decline of persistence (mostly 
due to the preceding dome flat) can be seen clearly. By the end of 
a visit the persistence declines to 50-70 counts per exposure. 

data cube and flagged in the mask image. These pixels 
are not used in any of the subsequent analysis. 

Pixels in the region of the Littrow ghostpq (^1.604 /rm, 
but it varies with fiber), are flagged in the mask but are 
not marked as bad in the image since the effect can be 
negligible depending on the spectrum. 


4.8. Persistence 

Teledyne H2RG arrays, like many IR detectors, ex¬ 
hibit persistence behavior in which a latent image of a 
previous exposure appears in subsequent image s but at a 


fraction o f the original source flux or stimulus (Smith et 
al. 2008a). In H2RGs, the persistence generally decays 


exponentially with time with the resulting persistence 
representing a small percentage of the stimulus, although 
it can sometimes take a long time to be released. 

Normal persistence is not significant for most 
APOGEE exposures. The total persistence accumulated 
(in a dark exposure) 1800s after a “stimulus” exposure is 
only ^60 counts even for stimulus well depths of 25,000 
counts. However, some regions of the APOGEE detec¬ 
tors — the top ^1/3 of the “blue” array and around the 
perimeter of the “green” detector — have unusually high 
levels of persistence (which we sometimes refer to as “su¬ 
perpersistence” ). In these regions the total accumulated 
persistence 1800s after an exposure is ^10-20% of the 
stimulus counts. 

APOGEE spectra can be contaminated in the high per¬ 
sistence regions in two ways: (1) by stars in the same 
fiber from preceding plate visits (most significant when 
the preceding star is bright), and (2) by calibration ex¬ 
posures, such as the “dome” flats taken before each plate 
visit for throughput calculations. While the persistence 
from the former have spectral features that will corrupt 
a stellar spectrum, the latter are generally featureless 


The Littrow ghost is formed by dispersed light being reflected 
from the detector arrays, reflectively recombine d by the VPH grat¬ 
ing. and Anally re imaged onto the detectors (|Burgh et al.||2007| 
[Wilson et al.|2015| >. 


Fig. 8 .— (a) The distribution of magnitude difference (AH) be¬ 
tween stars in the blue superpersistence region and the star in the 
preceding plate visit (in the same flber) where a positive differ¬ 
ence indicates that the preceding star was brighter, (b) The rela¬ 
tive blue excess of the stellar spectra, (blueo 6 s-blue^o(i)/t>lue^o(i 5 
where blue^o^Z is the expected (model) blue ffux and is the product 
of the observed green ffux and an empirical relation (found from the 
normal persistence region) of the blue/green ffux ratio as a function 
of J-Ks color (1.092—0.0737[J-i<rsl)- The relative blue excess sig- 
nifles the level of persistence relative to the stellar spectrum. The 
superpersistence region is shown in red and the normal persistence 
region in blue. The line indicates the median of all stars at the 
same magnitude difference and the error bars the robust standard 
deviation. 


and act as a veiling component. Eortunately, APOGEE 
uses a fiber management scheme in which fibers are des¬ 
ignated for “bright”, “medium” and “faint” stars (al¬ 
though each group still spans a significant magnitude 
range) to help reduce cross-contamination between spec¬ 
tra on the detectors (i.e., to minimize the “spatial” wings 
of a bright star cont aminating the spectru m of an adja¬ 
cent faint star; see [Majewski et al.j 2Q15| for more de¬ 
tails). This observing scheme hber management scheme 
also helps reduce the relative effect of persistence from 
preceding stars because the magnitude differences are 
smaller. 

Eigurej^ shows the median continuum in the blue de¬ 
tector for sky fibers (which should show persistence from 
previous exposures) versus exposure number for many 
science visits. This figure illustrates the level (and tem¬ 
poral behavior) of the persistence (for which the accu¬ 
mulated charge is dominated by the dome fiat immedi¬ 
ately preceding a science exposure), which mostly affects 
the fainter stars. While featureless (gray) persistence 
(as expected from the dome flat) will not significantly 
impact the measured radial velocities, it will affect the 
relative depth of the spectral absorption lines and there¬ 
fore the derived abundances. The effec ts of persistence 
on DR 12 abundances are described in IHoltzman et al.l 
(2015). Eigure shows the distribution of magnitude 
ditterence {AH=H 2 -Hi) between a star and the previ¬ 
ous star in the same fiber (from a different plate) while 
the “relative blue excess”, a measure of persistence rela¬ 
tive to the stellar spectrum, is shown in Eigurej^. Even 
when the preceding star is significantly fainter there is 
excess flux in the blue (^10%), likely due to persistence 
from the dome flat. As AH increases (i.e., the preceding 
star gets relatively brighter) the flux ratio also increases 
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Relative Blue Excess 


Fig. 9.— The distribution of relative blue excess for the super¬ 
persistence region (red) and a similarly sized portion of the normal 
persistence region (blue). The flux ratios are higher in the super¬ 
persistence region due to persistence from the dome flat and the 
preceding star in the same flber. The superpersistence distribution 
peaks at ~0.174 but has a long tail to higher values. 


due to the persistence from the preceding star. The stel¬ 
lar and dome flat persistence reach parity at AiJ^+0.10 
mag and a relative blue excess of 1 (where the persis¬ 
tence is as bright as the star itself) is reached at Ai7^+6 
mag. The median relative blue excess in the superper¬ 
sistence region is ^17.4% but with a long tail to higher 
values (Eigurej^ and, as might be expected, the faintest 
APOGEE stars are the most affected by the persistence 
(Eigurep^. 

Our initial investigation of the APOGEE superper¬ 
sistence behavior indicates that the persistence from a 
single stimulus (e.g., a flat exposure) is well de scribed 
by a double- exponential in time (as also seen by Smith 
et al. 2008b) with the timescales being fairly constant 
across pixels and stimulus (^120s and ^1700s). The ex¬ 
ponential amplitudes are a complex function of fluence 
(or count rate) as well as total exposure time, and the be¬ 
havior changes once saturation is reached. However, per¬ 
sistence becomes significantly more complex once multi¬ 
ple stimulus exposures are taken in a row. While this 
is the information that is required to properly correct 
APOGEE data for persistence (i.e., the persistence due 
to the multiple previous stimuli), we have not yet fully 
characterized this behavior. 

While there is evidence that persistence can be cali¬ 
brated out in some circumstances (e.g, the HST WEG3 
team has developed an algorithm that removes ^90% of 
the persistence in WEG3/IR image, the APOGEE 
reduction pipeline for DR12 contains no correction for 
persistence. This is largely due to the complexity of 
the problem. Our initial investigation has shown the 
WEC3 prescription to be inadequate for the APOGEE 
detectors. Work is proceeding to fully characterize the 
APOGEE persistence behavior and to develop an algo¬ 
rithm to correct for persistence in later data releases. Eor 
DR12, spectra in the persistence region and with jumps 
in the spectral continuum from green to blue are flagged 


http: //www. stsci. edu/hst/wf c3/ins_perf ormance/ 
persistence/ 


6 8 10 12 14 

H 



Fig. 10. — The distribution of relative blue excess with H- 
magnitude for the superpersistence region. The connected red dots 
show medians in 0.5 mag bins of H. The persistence is most sig- 
niflcant for the faintest APOGEE stars. 


(see |Holtzman et aLl|2015|) and should be treated with 
thT 


caution. Of course, the red and most of the green por¬ 
tions of these spectra are uncontaminated by superper¬ 
sistence and can be used for science. Looking forward, 
persistence is expected to be less of a concern in the 
SDSS-IV/APOGEE-2 survey because the blue detector 
was replaced in summer 2014 with a detector exhibiting 
normal persistence behavior (although the original green 
detector with higher persistence around the perimeter 
still remains). 


4.9. Output 

All of the above tasks are performed by the IDL rou¬ 
tine AP3DPROG program in the reduction pipeline. The 
output files are called ap2D-[abc]-ID8.fits (the abc denot¬ 
ing the red/green/blue detectors, respectively) and have 
three data extensions, one each for the flux, errors, and 
bitwise pixel mask, and each with a size of 2048x2048. 
These files are output to the spectro/v#/red/MJD5/ di¬ 
rectory. 


5. AP2D: EXTRACTION TO ID SPECTRA 

Once the 2D images are created, the next step is to ex¬ 
tract the 300 fiber spectra, correct them for fiber-to-fiber 
throughput variations, and wavelength calibrate them. 
This step is performed by the IDL AP2DPROG routine. 

5.1. Extraction 

The spectra are recorded roughly along rows on the 
three detectors. The standard extraction method uses 
an empirical measurement of the spatial profile of each 
trace at each column of each detector. Use of an em¬ 
pirical spatial point spread function (PSE) is possible 
because this PSE is quite stable over the course of an 
exposure sequence (visit), because the instrument is sta¬ 
tionary and fiber fed. Over a longer period, the traces 
do move by a fraction of a pixel as the weight of the LN 2 
tank suspended below the cold plate varies throughout 
a night due to LN2 depletion and then refill. However, 
this longer period drift is tracked using individual “dome 
flat” exposures that are taken at the end of every expo¬ 
sure sequence. These exposures are taken of the telescope 
mirror covers illuminated using a continuum source, and 
are used both for spatial profile construction as well as for 
mapping the system fiber-to-fiber throughput variations. 















10 


NIDEVER ET AL. 



Fig. 11.— An example of the PSF model (red) of a stellar spec¬ 
trum (black) along a section of one column (X=400) in the green 
detector. The residuals (data—model) are shown in the lower panel 
as well as ± the Poisson noise indicating that the model is satis¬ 
factory. 


The 300 spectra are separated by rQug:hly 6-7 p i xels on 
average. In detail, as described by Wilson et al. (2012), 
the spectra are grouped in 10 blocks ot 30 hbers each, 
with slightly larger gaps between these groups. 

To allow accurate measurements of the wings of the 
PSF, data are periodically taken using the so-called 
“sparse pack” calibration channel, in which only 50 of 
the 300 fibers are populated. These provide widely sep¬ 
arated spectra from which the PSF can be determined 
accurately. 

From the “dome” fiats, a rough correction for scattered 
light is made by subtracting the mean value measured at 
the top and bottom of the chips (specifically, rows 5-10 
and 2038-2042). The resulting image is then smoothed 
in the wavelength direction with a boxcar filter of width 
50 pixels, to increase S/N, given that the PSF is not ex¬ 
pected to change significantly on this scale (there is very 
little tilt of the traces against detector rows). At each 
column, the observed profile is then tabulated for each of 
the 300 traces. Since there is some overlap between adja¬ 
cent traces, the light at each pixel is distributed between 
the two surrounding traces using information from the 
sparse pack calibration frames. Based on the distances 
of a given pixel from each of the surrounding traces, the 
sparse pack calibration is used to determine the relative 
fraction of the two contributions by looking at values at 
these distances (on the appropriate side of the PSF) in 
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Fig. 12 .— An example of fiber-to-fiber throughput variation cal¬ 
ibration functions for all 300 fibers and three detectors. The rms 
variations are 9.7, 12.2, and 11.3% for the red, green and blue 
detectors respectively. 


Wavelen^^h^(|a.m) 



Fig. 13.— The response curve calibration spectrum. Each spec¬ 
trum is multiplied by this function to roughly calibrate the relative 
fluxes. 


the nearest trace in the tabulated sparse pack data. 

This method assumes that the profile is sufficiently lim¬ 
ited in spatial extent so that only adjacent traces con¬ 
tribute to the light at any given pixel. With this as¬ 
sumption, the brightness of any given trace is directly 
coupled only to the brightness of adjacent traces. Since 
the solutions of neighboring triplets of fibers are linked 
together, the solution of all fiber fluxes becomes a linear 
algebra problem with a tridiagonal matrix, solving for 
300 individual fluxes at each column. This problem is 
solved b y using the trid iagonal matrix or “Thomas” al¬ 
gorithm (Thomas||1949). Errors are propagated through 
the extraction. An example of a model PSF fit to the 
APOGEE data is shown in Eigure IE _ 

The AP2D code also includes options for three other 
extraction methods: boxcar extraction, and simultane¬ 
ous fitting using both a simple Gaussian and a more 
complex functional (but still parametric) form for the 
PSE. Euture work might include developing these fur¬ 
ther. The most significant future development, however, 
would probably be to go to a full two-dimensional ex¬ 
traction that accommodates a 2D PSE that is not sepa- 
rable in rows and colnni ns (i.e., “spectro-perfectionism”, 
Bolton fc Schlegel]|2QlQ ). 


5.2. Fiber-to-fiber throughput variation and response 
eurve eorreetion 

The “dome” fiats for each plate are then used to correct 
the 300 extracted spectra for moderately wavelength- 
dependent fiber-to-fiber throughput variations. These 
are measured on an individual plate basis because of the 
possibility that the relative throughputs vary depending 
on the details of each specific coupling of the long fibers 
from the instrument to the short fibers in the plug plates 
via the gang connector. Eigure [T^ shows an example of 
a flux calibration file. The fiber-to-fiber variations are at 
the ^10% level. 

Accurate fiber-to-fiber throughput correction is desired 
to enable the possibility of accurate sky subtraction us¬ 
ing separate sky fibers. The correction reduces the rms 
variation of the flux in sky fiber lines across the plate 
from ^12% to ^5% and only ^1% around a smoothly 
varying 2D spatial polynomial fit to the flux variations. 

After the fiber-to-fiber throughput calibration, each 
spectrum is corrected by a wavelength-dependent spec- 
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(normal science exposures, sky flats, and Any Star Down 
Any Eiber [ASDAPj^ exposures) the night sky airglow 
emission lines are used for this correction. A zero-point 
offset in the pixel positions (not the wavelength) is deter¬ 
mined for each fiber separately. Eirst, Gaussians are fit 
to all the bright emission lines in each spectrum. Then, 
a first guess zero-point offset is determined by cross- 
correlating a model airglow spectrum (using Gaussians 
with heights of unity and a standard APOGEE wave¬ 
length solution) with a model spectrum of the measured 
lines (also using Gaussians with heights of unity). Next, 
the zero-point estimate is used to match the measured 
lines with known airglow lines. A new wavelength solu¬ 
tion is determined allowing only the pixel zero-point to 
vary. Einally, a line is fit to all the fiber zero-point shifts 
versus fiber number and this fit is used for the zero-point 
offsets in the wavelength solution. 


Fig. 14. — Residuals of Thorium-Argon-Neon and Uranium-Neon 
line fits around the wavelength solution versus position along the 
three detectors (color-coded by position in the spatial dimension). 
The global rms is 0 . 028 A and dominated by the measurement un¬ 
certainty in the line centers. Connected black crosses show medians 
in bins of 200 pixels. Systematics in the binned residuals are at 
the '^0.01-0.02A level. 

tral response function to apply an approximate relative 
flux calibration. The relative response curve was created 
using the spectrum of a blackbody of well-known temper¬ 
ature (110 °C and 150 °C) and is used for all APOGEE 
spectra (Eigurep^. 


5.3. Wavelength Calibration 

The correspondence between wavelength and pixel is 
determined using exposures of ThArNe and UNe hollow- 
cathode lamps. The wavelength solution varies from 
fiber-to-fiber because of the optical design and also be¬ 
cause of the exact placement of each fiber in the pseudo¬ 
slit. The former causes variations over large scales, while 
the latter causes a flber-to-flber shift. 

To derive the wavelength solution, Gaussians are fit 
to all detectable lines (4cr above the background) in all 
calibration spectra. Next, these Gaussians are matched 
up to known lines. The information for all the lines in 
the 300 fibers, three arrays, and multiple exposures (nor¬ 
mally one ThArNe and one UrNe) are combined into one 
list. A 5th order polynomial (A versus pixel) and two de¬ 
tector gaps are fit to each fiber separately. A robust 
linear fit (to allow for the possibility of small rotations 
of the detectors relative to each other) is made to the 
derived detector gaps as a function of fiber number to 
obtain more accurate values. Einally, the 5th order fits 
are redone holding the chip gaps fixed at the fitted val¬ 
ues. The residuals to these fits are on the order of ^0.03- 
0.04A (^0.1 of a pixel) and dominated by the Gaussian 
line-fitting errors. Systematics in the residuals are at the 
^0.01-0.02A level (Eigure 14). The temporal variations 
over year timescales is ^0.015A, which demonstrates that 
the instrument is very stable. 

A wavelength solution from an apWave calibration file 
is applied to each extracted spectrum. This wavelength 
solution still needs to be corrected for slight differences 
between the science exposures and the calibration frames 
due to different dither po sitio ns and potential changes in 
the optics over time (see ^6.2). Eor “on-sky” observations 


5.4. Output 

The output files for AP2D are called aplD-[abc]- 
IDS.fits and have four extensions containing the flux, er¬ 
rors, bitwise pixel mask, and wavelength array (in A), 
each with a size of 2048x300. A model of the 2D image 
is also output with names of ap2Dmodel-[abc]-ID8.fits. 
These files are output to the spectro/v#/red/MJD5/ di¬ 
rectory. 

6 . APIDVISIT: VISIT STAGE 

After extraction, the spectra are sky and telluric cor¬ 
rected and then the separate exposures (at different 
dither positions) are combined into one well-sampled 
spectrum per fiber. This is performed in the APIDVISIT 
stage with the APIDVISIT routine. 

As explained in Section each plate “visit” generally 
consists of eight ^500 sec (47-read) exposures taken as 
two ABBA dither position sequences (A and B being 
roughly 0.5 pixels apart in the spectral dimension). 

Sky subtraction and telluric absorption correction are 
done on an exposure-by-exposure basis because of the 
possibility of sky variation from one exposure to another, 
since the sky is known to vary on short time scales. This 
leads to some complications because of the mild under¬ 
sampling of the spectra on the blue end, preventing sim¬ 
ple resampling of sky and telluric fibers to yield a cor¬ 
rection for each object fiber. 

6.1. Dither Shift Measurement 

Due to slight undersampling at the bluer wavelengths 
the APOGEE exposures are taken at two different spec¬ 
tral dither positions shifted by roughly 0.5 pixels from 
each other. The actual positions at which the exposures 
are taken are not known precisely enough a priori for 
accurate combination of spectra taken at two different 
dither positions. Therefore, the dither position of each 
exposure is measured from the actual data. This is done 
relative to the first (reference) exposure in two ways: 
(1) cross-correlation of spectra, and (2) shifts of airglow 
emission lines. Eor the cross-correlation, the (non-sky) 
spectra are first normalized (using a 100-pixel median fil¬ 
ter), then cross-correlated against the first exposure spec¬ 
trum, and finally a Gaussian is fit to the cross-correlation 

ASDAF exposures are special exposures of bright stars (often 
“standard” or calibration stars) placed on a specific fiber. 
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Fig. 15.— A diagram illustrating the LSF variations across the 
detectors. Each profile shows ±7 pixels of the model LSF convolved 
with the pixel size. 

peak to find the best shift. This is done separately fiber- 
by-fiber and array-by-array. Finally, a robust mean is 
calculated of all individual 900 shifts measured. In the 
emission line method, all bright emission lines in all fibers 
are fit with Gaussians, then matched with their corre¬ 
sponding lines in the first exposure, shifts calculated, and 
a robust mean measured. The pipeline currently uses 
only the cross-correlation technique which gives average 
formal errors of ^0.005 pixels. The measured dither 
shifts are written to the header for later use in the dither 
combination process. 

6.2. Line Spread Profile (LSF) 

To accomplish sky correction (emission and absorp¬ 
tion) on undersampled data, as well as for spectra in 
which the LSF varies across the detector, the sky cor¬ 
rections are forward modeled using constraints from the 
observed sky and telluric fibers. This requires an accu¬ 
rate measurement of the LSF as a function of wavelength 
and fiber. 

The line spread function (LSF) is modeled as a sum of 
Gauss-Hermite function^^ (which form an orthonormal 
basis) and a wide-Gaussian for the wings. The convolu¬ 
tion of the model LSF with the pixel size is taken into 
account. LSF parameters are determined using airglow 
lines in the nightly skyflats. Each fiber and array are 
fit separately. The LSF parameters are allowed to vary 
slowly with wavelength by means of a low-order polyno¬ 
mial as a function of pixel. The model LSF is always 
normalized to an area of one. Figure illustrates the 
typical variations of the LSF across the three detectors 
and the fact that the LSF is slightly undersampled in the 
blue. 

As a demonstration of the LSF stability. Figure pT] 
shows the temporal variation of the Gaussian FWHM 

We use the probabilists’ polynomials, which are given 

by: Hen{x) = ( —l)^e^ ’ https://en.wikipedia.org/ 

wiki/Hermite_polynomials 


of one Thorium-Argon-Neon line per detector (and five 
fibers). While there are some small-scale variations, the 
rms of the FWHM values per fiber are ^1-2%, which in¬ 
dicates that the instrument has been remarkably stable 
over three years. 

The traditional measurement of resolving power 
(R=A/AA; using a direct-measured FWHM for A A) 
gives an average of R=22,500 for APOGEE spectra, but 
there are ^10-20% variations across the detectors (spec¬ 
trally and spatially) as seen in Eigure 17 However, the 
APOGEE LSE is non-Gaussian and, tnerefore, a tradi¬ 
tional “R”-value is not necessarily a good description of 
apogee’s resolution. 

6.3. Sky Subtraetion 

The observed spectra are contaminated by night sky 
emission airglow lines (mostly OH) from the Earth’s at¬ 
mosphere, and sky continuum which is most often dom¬ 
inated by moonlight (reflected sunlight) and greatly en¬ 
hanced on cloudy nights. We have found that light pol¬ 
lution from nearby El Paso (Texas) is not a significant 
component of our sky spectra. 

The sky spectrum can vary spatially across our 3° EOV 
and temporally during our hour “visit” of the held. 
Therefore, each plate has ^35 fibers designated for 
“blank” sky that can be used to subtract the sky spec¬ 
trum from the science fibers. The sky fiber positions are 
chosen in ^17 spatial zones to giv e them a fairly uniform 
distribution across the plate; see Zasowski et al. (2013) 
for details. 

The current pipeline is using a temporary, and sub- 
optimal, sky subtraction method. Eor each science fiber 
the four nearest (on the sky) sky fibers are found. Eirst, 
the spectra from these four sky fibers are resampled onto 
the wavelength solution of the science fiber (using cubic 
spline interpolation). Next, each of the four sky spectra 
are cross-correlated with a continuum subtracted (150 
pixel median Alter) science fiber spectrum and shifted 
accordingly to fix any errors in the wavelength solutions. 
An emission line scaling factor is then calculated for each 
sky spectrum relative to the science spectrum airglow 
lines using pixels with sky emission greater than 10 x the 
noise. A weighted average sky spectrum is then created 
from the scaled sky spectra (with sky continuum added 
back in) using l/(sky distance)^ for the weighting, and 
subtracted from the science spectrum. Eigure \T8\ shows 
an example of a single exposure sky fiber spectrum af¬ 
ter sky subtraction (the absolute value; the sky fiber in 
question was not used to determine the sky spectrum, 
only the neighboring sky fibers were) and the pipeline 
noise model for comparison. Even this sub-optimal sky 
subtraction works fairly well with the large majority of 
the residuals below three times the noise level. 

Possible improved methods include modeling of the air¬ 
glow lines with 2D spatial polynomial fitting of airglow 
family fluxes, and principal component analysis (PGA). 
These methods are currently under investigation and de¬ 
velopment. 

It should be noted that the airglow lines are so bright 
that little scientific benefit can be gained from the stel¬ 
lar spectrum “underneath” them, since even with perfect 
subtraction the Poisson noise from airglow will dominate 
over (wash out) any stellar spectrum. The far wings of 
the airglow lines should be salvageable, but the most crit- 
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Fig. 16. — A demonstration of the temporal stability of the LSF. The Gaussian FWHM (in pixels) of one ThArNe line per detector and 
five fibers versus MJD. The FWHM values are stable at the ~l-2% level over three years. 


ical portion of sky subtraction is the removal of the sky 
continuum, which would otherwise distort line depths in 
normalized spectra and their subsequent analysis. Poorly 
subtracted airglow lines can also adversely affect the RV 
determination. 


6.4. Telluric Correction 

The NIR i^-band hosts a number of atmospheric (tel¬ 
luric) absorption lines (from H 2 O, CO 2 , and CH 4 ; see 
Eigure 19) that contaminate a significant fraction of the 
APOGEE spectral range (^20%) and need to be cor¬ 
rected. As with the airglow spectrum the telluric ab¬ 
sorption spectrum may vary spatially and temporally for 
our observations. Each plate has ^35 fibers designated 
for “telluric” (hot) stars that can be used to ascertain the 
telluric absorption for the science fibers. A procedure is 
used to pick telluric fibers in spatial zone s similar to that 
used for sky fibers (Zasowski et al. 

A three-step process is used to correct for telluric ab¬ 
sorption: ( 1 ) telluric absorption model fitting to the hot 
star spectra, (2) 2D polynomial spatial fitting of the tel¬ 
luric species scaling parameters across the plate, and (3) 
construction of the model telluric absorption spectrum 
for each science fiber using the model telluric spectra, 
2 D polynomial fitting parameters, and the known LSE 
of the science fiber. 

The L BLRTIVg model atmosphere code ( [Glough 


aLl|2005|) was used to create a grid of high-resolution 
model telluric species spectra for GO 2 , H 2 O and GH 4 , 
individually, using the US Standard 1976 atmosphere for 
the altitude of APO. Eor GO 2 and GH 4 , four different 
scale factors that parameterize the strength of the fea¬ 
tures were used (0.5, 1.0, 1.5, and 2.0 cm); for H 2 O, four 
different precipitable water columns (0.75, 1.5, 2.25, and 


http://rtweb.aer.com/lblrtm.html 


3.0) were used. Spectra were calculated at seven different 
airmasses, from 1.0 to 2.5, spaced by 0.25. 

Eor each plate, we interpolated in airmass within this 
grid to get four model spectra for each absorption species. 
Eor each telluric star on a plate, we find the scale factor 
of each model that best matches each spectrum, and de¬ 
termine which of the four scaled models for each species 
provides the best fit. We then adopt the model for each 
species that best fits the majority of the stars, and re¬ 
fit each telluric spectrum with the same model to get a 
self-consistent set of scale factors across the field. 

The fits are performed by adopting the scale factor 
that yields the minimum RMS in all pixels where the 
telluric lines for that species are dominant (pixels must 
be within 5 pixels of a telluric species line with strength 
greater than 1 % in the convolved model telluric spec¬ 
trum). To obtain a reliable telluric scale factor the stel¬ 
lar continuum (including the wide hydrogen absorption 
lines) needs to be removed. Therefore, an iterative pro¬ 
cess is used. The stellar continuum is calculated using 
a median filter (100 pixels wide) of the stellar spectrum 
corrected with the current best-fit model telluric spec¬ 
trum (after the first iteration). Stellar continuum and 
scale factors are calculated until the solution converges. 
Since the line spread profile (LSE) varies from fiber to 
fiber, the high-resolution model species spectra are con¬ 
volved with a separate LSE (from the apLSE calibration 
file) for each fiber. 

At this point the species scalings for each hot star are 
determined, and a 2D spatial polynomial model is fit to 
the ^35 scalings for each species separately. Eor GO 2 
and GH 4 we use a linear model; since these species are 
thought to be well-mixed in the atmosphere, the linear 
trend is included to account for variations in airmass 
across the field. Eor H 2 O, which could have significant 
spatial structure, we use a quadratic surface. The RMS 
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Fig. 17.— Maps showing how the resolving power (R=A/AA) varies across the three detectors. 


around the fit is normally ^1-2%. Eigure 20 shows an 
example of a 2D polynomial fit to a species scaling with 
significant spatial variations (most exposures show much 
less variations). 

Einally, the model telluric spectra are calculated for 
each science spectrum, using the 2D polynomial fit coeffi¬ 
cients to obtain the three species scalings for the position 
of the science object. Then, the high-resolution model 
telluric species spectra are convolved with the LSE of the 
science fiber and scaled appropriately with the species 
scaling for that fiber. The science spectrum is divided 
by the final convolved telluric spectrum. An error in the 
telluric model for each object is computed by taking the 
RMS scatter in the 2D polynomial fit of each species scal¬ 
ing and propagating the errors forward into the model 
telluric spectrum. The average error in the telluric cor¬ 
rection is roughly ^1-2 % of the stellar continuum. It 
is worth noting that this procedure is astrophysically in¬ 
correct as the telluric absorption occurs before the light 
goes through the spectrograph and gets convolved with 
the instrumental LSE. However, this often-used approx¬ 
imation works well enough (^1%) for our needs. 

Eigureshows an example of an APOGEE spectrum 
of a hot star before and after telluric correction. The 
residuals after removing the broad stellar absorption fea¬ 
tures are small at 0.0081 or less than a percent. Some sys- 
tematics in the residuals, such as around the CO 2 band 
on the red side of the blue chip, are still present. Since 
the LSE varies only slowly with time it is not necessary 
to convolve the high-resolution telluric model spectrum 
for each exposure separately. Therefore, to save time in 
the entire telluric correction step, the three original high- 
resolution telluric model spectra are pre-convolved with 
the LSE of all 300 fibers and saved as the apTelluric cali¬ 
bration product. Currently, only a single LSE calibration 
file is used for all the APOGEE reductions, because the 
LSE shows very little temporal variation (see Eigure 

After sky subtraction and telluric correction the ‘cor¬ 
rected” frame is written to disk. The apCframe-[abc]- 
IDS.fits files have extension of 1-flux (electrons), 2-error 


(electrons), 3-bitwise flag mask, 4-wavelength (A), 5-sky 
(electrons), 6-sky error (electrons), 7-telluric absorption, 
8-telluric error, 9-wavelength coeficients, 10-LSE coeffi¬ 
cients, and 11-plugmap structure (binary table). The ar¬ 
rays in extensions 1-8 have sizes of 2048x300, the wave¬ 
length coefficients 300x14, and LSE coefficients 300x27. 

6.5. Dither Combination 


As mentioned in §6.1 
APOGEE exposures 


are 


due to undersampling the 
taken at two different dither 
positions that must be combined to create well-sampled 
spectra. This is performed in two steps: (1) the ex¬ 
posures (a total of Nexp) are paired up and interlaced 
to create well-sampled spectra, and (2) the Nexp/2 well- 
sampled spectra are co-added to create one “visit” spec¬ 
trum per object. 

The exposures are paired up, one dither position (A 
and B) per pair to create a series of approximately equal¬ 
s'/ A" pairs. Each exposure spectrum is separately nor- 
malized/scaled using a median filtered (width of 501 pix¬ 
els) version of the spectrum (each array separately) to 
remove any variations in net flux due to differences in 
seeing or throughput between the two exposures (only for 
objects, not sky). Using the previously measured dither 
shifts ( §6.1| ), the two scaled spectra are then combined 


with the smc-interlace equations from Bracewell 

^'native” sea 


(1999) 
e (for 


onto a pixel scale twice as fine as the 
the single exposures). All pairs share the same final pixel 
scale so that the spectra are only resampled once. The 
flux level of the well-sampled spectrum is then rescaled 
to the average of the scaling arrays of the two individual 
exposures. 

In the second step, the Nexp/2 well-sampled spectra 
(normally 4 per visit) are co-added. The flux level of the 
spectra are again scaled using a median filter. The scaled 
spectra are combined using a weighted mean (with 5 a 
outlier rejection) where the weights are computed either 
on a spectrum-to-spectrum or pixel-by-pixel basis. The 
pipeline currently uses the pixel-by-pixel weighting. The 
final combined spectrum is then rescaled using the sum 
of the scaling arrays of the individual spectra. 
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Fig. 18. — Example of a single exposure APOGEE sky fiber spectrum (exposure 02870014, fiber 258) after sky subtraction (the absolute 
value; the sky fiber in question was not used to determine the sky spectrum, only the neighboring sky fibers were) with the pipeline noise 
model for comparison. Most sky line residuals are below 3x the noise model. 


6.6. Absolute Flux Calibration 

The final combined object spectra are roughly flux cal¬ 
ibrated (on an absolute flux scale) using the 2MASS H- 
band magnitudes. Eor each star we separately deter¬ 
mine an absolute flux correction factor, Cabs^ to allow for 
differences in airmass, fiber-centering errors, and seeing 
variations. This multiplicative factor is used to convert 
the observed spectrum in electrons to calibrated physical 
units. 


Ea = CabsSi (1) 

where S is the observed spectral flux in e“ at pixel i. 
The correction factor can be calculated with. 


Cabs = MEDIAN {Si) 


( 2 ) 


where Fa,o is the 2MASS zero-point for isophotal 
monochromatic light (1.33xlQ ~^^ W/cm^//im) fo r a 0th 
magnitude star in the i7-band (Cohen et al. |2003P j Af- 
ter applying the absolute flux calibration the spectra are 
in physical units of ergs/s/cm^/A. 


6.7. Output 

After the absolute flux calibration step the apPlate- 
[abc]-PLATE4-MJD5.fits and apVisit-PLATE4-MJD5- 
EIBERID3.fits files are written to disk. The apPlate 
files contain all 300 spectra while the apVisit files (^265 
of them) are for single object spectra (no sky spec¬ 
tra). The fluxes and errors are stored in units of 10“^^ 
ergs/s/cm^/A. 

http://www.ipac.caltech.edu/2mass/releases/allsky/ 
doc/sec6_4a.html 


Note, the EiberlD is not the same as the IDL index 
in the data arrays. EiberlD=300 is the first spectrum in 
the data arrays (bottom) and EiberID=l is the last one 
(top). The conversion is Eiberlndex=300—EiberlD. 

Many of the important object parameters includ¬ 
ing name, coordinates, 2MASS magnitudes, APOGEE 
targeting flags, date observed, apVisit filename, ra¬ 
dial velocity, best-fitting template parameters, and 
more are saved in a EITS binary table called 
apVisitSum-PLATE4-MJD5.fits. in the field directory 
spectro/v#/f ields/L0C4/. These are useful for quickly 
accessing the output parameters of the pipeline. 

7. APSTAR: OBJECT STAGE 

Most of the stars are observed in several different visits 
to enable detection of radial velocity variation to identify 
binaries and to accumulate the necessary S/N. After 
multiple visits, a combined spectrum is made to provide 
the highest S/N individual spectra for each star. 

The output combined spectra for all APOGEE objects 
are placed on the same rest wavelength pixel scale, with 
a constant dispersion in log A, using 

\ogXi = 4.179-1-6 X 10-^i 

with 8575 total pixels (i = 0 to 8574), giving a rest wave¬ 
length range of 15100.8 to 16999.8 A. The dispersion was 
chosen to provide approximately 3 pixels per resolution 
element, although the resolution varies over the full wave¬ 
length range. 

To do the combination, each rest frame (i.e., the 
Doppler shift is removed) visit spectrum is sampled on 
this final wavelength scale using sine interpolation. Since 
the visit spectra are all dither combined, they are well 
sampled over the entire range; in fact, they are signif¬ 
icantly oversampled at the long wavelength end. The 
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Telluric Absorption Spectra 



Fig. 19.— The atmospheric telluric absorption features in the APOGEE spectral window due to CH 4 (green), CO2 (red) and H2O (blue). 


Species H20 Species H20 



0.410 0.422 0.435 0.447 0.459 0.471 0.484 



Fig. 20.— Example of a 2D polynomial fit to the telluric species 
scaling for H2O (exposure 03770010). C 'H ^'^e coordinates 
relative to the center of the plate in right ascension and declination, 
respectively. This is a fairly extreme case to show how the 2D 
fitting handles large spatial variations. Most cases are flatter. 


sine interpolation takes this into account by using a chip- 
dependent EWHM, conservatively adopted to be 5, 4.25, 
and 3.5 (dithered) pixels in the red, green and blue chips, 
respectively. This effectively filters out noise at higher 
spatial frequencies. 

After a rough continuum normalization (using wide 
boxcar and median filters), the resampled spectra are 
then combined using a weighted mean, with weights 
calculated on both a pixel-by-pixel and a spectrum-by¬ 
spectrum basis (where the weight is the inverse square 
of the normalized error). The resulting weighted average 
spectrum is then multiplied by the average continuum of 
the individual visit spectra. The determination of radial 
velocities and the spectral combination are done itera¬ 


tively and is described in more detail in Section 8.2.3[ 

Eigure 22 shows the S/N values (per half resolution 
element) using the pipeline noise model and the variance 
from multiple visits (for 9548 stars with 6 visits). The 
empirical values are systematically low compared to the 
noise model estimates in the high S/N regime (i.e., bright 
stars), which indicates that we are limited by systematics 
(at the ^0.5 % level) for these stars. However, we can 
easily achieve the S/N^lOO required for the survey. 

A combined LSE is created by taking a weighted av¬ 
erage (using the same weights as above) of the indi¬ 
vidual visit LSE arrays on the final apStar wavelength 
scale. The model LSE (a sum o f spa tially-varying Gauss- 
Hermite functions; see Section 6.2) is fitted to the em¬ 
pirical 2D LSE array and the coemcients of this approx¬ 
imation are also saved. 


7.1. Output 

There is an output apStar-2MASSID.fits file for each 
unique object that contains the two combined spec¬ 
tra (pixel-by-pixel and visit-by-visit weighting) as well 
as the individual visit spectra resampled on the ap¬ 
Star wavelength scale. The final, weighted LSE and 
the Gauss-Hermite coefficients are saved in apStarLSE- 
2MASSID.fits files. Summary information for all stars in 
a given field are saved in the apEield-LOG4.fits files. 

8. RADIAL VELOCITY DETERMINATION 

APOGEE radial velocities (RV) are derived at both 
the Visit and the Object stages. The main steps are: 

• As each visit is reduced, an RV estimate is de¬ 
termined by cross-correlating the visit spectrum 
against a grid of synthetic spectra. This provides 
an “estimated RV” for the visit, which is stored in 
the apVisit files, but not subsequently used since it 
was found to be unreliable sometimes. 

• “Refined” radial velocities for each visit are derived 
when the visit spectra are combined in the object 
stage. This is done in three steps: 

1. Relative radial velocities are determined using 
the combined spectrum as the spectral tem¬ 
plate. This is done iteratively because the rel- 
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Fig. 21.— Example of an APOGEE spectrum of a hot star before and after telluric correction. (Top) Original spectrum in counts with 
telluric corrected spectrum below and the median-filtered stellar “continuum” (red). (Bottom) The “residuals” (red; telluric corrected 
spectrum divided by the stellar “continuum” minus 1) with an RMS of 0.0081. 


S/N in combined spectrum (half res element, 6 visits) 



Fig. 22.— The S/N per half resolution element of pixels not af¬ 
fected by telluric absorption or sky lines derived using the pipeline 
noise model (black) and the variance in multiple measurements 
(red) for 9548 stars with 6 visits. The lines show medians in bins 
of 1 mag. The empirical S/N values are consistently low compared 
to the noise model values for high S/Ns, which suggests that we 
are limited by systematic errors at the ~0.5% level for these stars. 

ative RVs and the combined spectrum depend 
on each other. 

2. The absolute radial velocity of the combined 
spectrum is determined by cross-correlating it 
with a grid of synthetic spectra spanning a 
large range of stellar parameters. 

3. The relative radial velocities for each visit and 
the absolute velocity of the combined spec¬ 
trum are then combined to produce absolute 
velocities for all visit spectra. 


lined spectroscopic binaries). This allows us to create 
a high-quality combined spectrum without even know¬ 
ing what type of object we are dealing with. However, 
the absolute RV is a critical science product and the fi¬ 
nal combined spectrum must be on the rest wavelength 
scale so that it can be properly compared to the large 
grid of synthetic spectra in the abundance pipeline (AS- 
PCAP). Therefore, the second step in the “refined” RV 
determination is to derive the absolute radial velocity of 
the combined spectrum against a small grid of synthetic 
spectra (the “RV mini-grid”). 

The various steps in the latter process are described in 
more detail below. 


8.1. Visit Stage Radial Velocities 

At the end of the ID Visit stage RVs are derived for 
all object spectra on the plate using the well-sampled 
apVisit files. This gives a first estimate for the RV (and 
therefore called “estimated” or EST RVs), but currently 
these values are not used later on in the processing. 
These were the only RVs avai lable pre-DRlQ and were 
used for science papers such as 
decided to keep this portion o: 


Nidever et al. (2012). We 
the code in the pipeline 
so we could compare RVs between the various methods. 
Eor completeness we briefly describe how these estimates 
RVs are derived. 

RVs are determined using a set of 96 synthetic template 
spectra (first generation “RV mini-grid”) that sparsely 
cover a large range in stellar parameters: 


• 3,500 < Teff < 25,000 K 

• 2.0 < log^ < 5.0 

• -2.0 < [Ee/H] < 0.3 


The latter scheme was employed because RVs derived 
from the combined spectrum (of the star itself) should 
be more precise than RVs derived from a small set of syn¬ 
thetic spectra (although there can be issues for double- 


The ASSeT ( [Koesterke et ah 2008 Koesterke 2009 


2012) spect rar synthesis package and Kurucz model at- 


mospheres (Castelli & Kurucz 2004 Meszaros et al. 2012) 
were used for cool and warni stars {/I'eff < 10,000 k) 




















































































18 


1.4 

1.2 

1.0 

0.8 

0.6 

0.4 

0.2 

0.0 


NIDEVER ET AL. 

apVisit-r5-4829-55733-028 2M20161309+3751193 



S/N=457.4 H=7.58 J-K3=0.66 ' ' ' '.'.-j 

- 

Observed 

Model 

_,_1_,_ 


1.65 

1.66 1.67 1.68 1.69 

i _^^_1_,_ ^^^^_,_. 

V,3i=-29'.82 VHeiio=-16.50+/-0.003 km/s.'.' 1 

1.59 1.60 1.61 1.62 1.63 1.64 

j ' ' ' ' 1 

=.1 

... 

.1 

Best Model Spectrum: [Fe/H]=0.00 TeJf=4,500 logg=3.50 x^=19.81 ' ' = 

' 1 ' ' 1 " 1 “ 

.1.1.1.1.1 = 


1.0 


1.0 

0.8 


0.0 


1.52 


1.53 


1.54 


1.55 

Wavelength (|im) 


1.56 


1.57 


1.58 


and Synspec dHubeny et al.||19 ^ iHubeny fc Lanz|2Qll|) 
TLUSTY model atmospher^ (Hifeny||l988['Hubeny ^ 


LanSlfl^ 




I 


Lanz & Hubeny 2UU3| 2(j(j7p 
UUU Kj. 'rhree steps are useo 


lor hot stars 
to derive the 


Fig. 23. — Example APOGEE normalized spectrum (black) with best-fitting template spectrum (dashed red) overplotted. This is the 
“estimated” RV-fitting done at the individual visit level. 

template is found using cross-correlation, a second RV- 
determination technique is used. The observed spectrum 
is split up into 45 pieces (^273 pixels each) and a sepa¬ 
rate RV derived for each piece using -minimization and 
the chosen template spectrum. In this simple forward- 
modeling technique, the only floating parameter is the 
Doppler shift. At a given Doppler shift the template 
spectrum is resampled onto the wavelength scale of the 
observed spectrum. One advantage of the “pieces” tech¬ 
nique is that an (internal) RV uncertainty can be cal¬ 
culated directly from the multiple RV measurements of 
the pieces. This technique works quite well for spectra 
of cool or metal-rich stars that have many narrow lines, 
but less well for hot stars. 

The final RV is chosen based on the calculated RV 
uncertainty of the two techniques. A final is calculated 
using the final, adopted Doppler shift and best-f itting 
tempiate. The barycentric correction (see Section 8.2.6 
beiow) is appiied to convert the RVs to the soiar system 
barycenter. An exampie APOGEE spectrum with the 
best-fitting tempiate is shown in Eigurel^ 


RVs: (1) continuum normaiization of observed and tem¬ 
piate spectra, (2) cross-correiation of observed spectra 
against aii tempiate spectra, and (3) y^-minimization us¬ 
ing the best-fitting tempiate from step 2. 

Normalization: Continuum normaiization of the 
spectra is criticai to obtain accurate RVs and especiaiiy 
important for hot stars with very few, but wide, spectrai 
features (i.e., Brackett hues). The spectrum for each of 
the three arrays is normaiized separateiy. Eirst, bad pix- 
eis and pixeis with bright airgiow iines are masked. Next, 
the 95^^ percentiie is caicuiated in 40 spectrai chunks (of 
^102 pixeis each), and then fit with a robust cubic poiy- 
nomiai. The spectrum is normaiized with this first esti¬ 
mate of the continuum, and the binning and poiynomiai 
fitting are repeated to remove some residuai structure. 
The finai continuum is the product of the two poiyno¬ 
miai fits. The same procedure is used to normaiize aii of 
the synthetic tempiate spectra. 

Cross-correlation: Before the normalized observed 
spectrum can be cross-correlated with the template spec¬ 
tra it is resampled (using cubic spline interpolation) onto 
the logarithmic wavelength scale of the template spec¬ 
tra. The resampled observed spectrum is then cross- 
correlated with the template spectrum, the best shift is 
found using the peak of the cross-correlation function, 
and y^ computed after shifting the template. This pro¬ 
cedure is repeated for all template spectra, and the best¬ 
fitting template is chosen based on the lowest y^ value. 
To refine the radial velocity determination, a Gaussian 
plus linear fit is performed on the peak of the cross¬ 
correlation function of the best-fitting template. 

y^-minimization method: After the best-fitting 


8.2. Object Stage Radial Velocities 

In the Object stage radial velocities are determined for 
all visits of a star together using a common RV template. 
The measurement of RVs and the spectral combination 
are performed iteratively as described below. 

8.2.1. Preparing the Speetra 
The spectra are “prepared” for cross-correlation by: 

• Pixel masking: Pixels marked as “bad” in the mask 
array or have sky lines in the sky array are masked 
out for the rest of the RV determination. 

• Continuum normalization: Each of the three chip 
spectra are normalized separately. The chip spec- 
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trum is separated into 40 chunks (covering approxi¬ 
mately 14 Angstroms each) and the 95^^ percentile 
pixel value is calculated for each chunk. A robust 
third-order polynomial is then fit to the chunk 95th 
percentile values. Einally, the spectrum is normal¬ 
ized (divided) by the polynomial fit. This is very 
similar to the Visit stage normalization method 
mentioned above. 

This preparation process is performed on both the ob¬ 
served visit spectra as well as the RV template spectra 
(observed combined or synthetic spectrum). 

8.2.2. Cross-Correlation 

All radial velocities are determined by cross-correlating 
a spectrum against a template spectrum. The spectra are 
on the same logarithmic wavelength scale (see Section]^ 
so that a Doppler shift is identical to a constant shift m 
the x-dimension. The spectra are “prepared” for cross¬ 
correlation by continuum normalization. A Gaussian is 
fit to the peak of the cross-correlation function to more 
accurately determine the best spectral shift. Einally, the 
shift and its uncertainty are converted to velocity units. 

8.2.3. Relative Radial Veloeities 

The relative radial velocities are determined by using 
the combined spectrum as the RV template. This is done 
iteratively, first determining the relative RVs and then 
creating the combined spectrum using the relative RVs to 
shift the visit spectra to a common (mean) velocity wave¬ 
length scale. Eor the first iteration, when no combined 
spectrum exists yet, the highest S/N visit spectrum is 
used as the template. Eor all subsequent iterations the 
combined spectrum is used as the template. Each itera¬ 
tion finds small shifts of the shifted and resampled visit 
spectra compared to the combined spectrum until the 
values converge. 

8.2.4. Absolute Radial Veloeities 

The combined spectrum after the relative RV step still 
has the mean RV of the star, which must be removed. 
The combined spectrum is cross-correlated against each 
synthetic spectrum in the RV mini-grid. Eor each syn¬ 
thetic spectrum the best RV and (of the shifted spec¬ 
trum) are derived. The spectrum with the lowest is 
chosen as the best-fitting spectrum and it’s RV is used 
as the absolute RV of the combined spectrum. Once the 
mean velocity is determined the visit spectra are com¬ 
bined one last time with the mean velocity removed so 
that the final combined spectrum is on the rest wave¬ 
length scale. 

The second generation RV mini-grid is composed of 
538 synthetic spectra that span a large range of stellar 
parameters: 

• 2,700 <Teff < 30,000 K 

• 0.0 < log^ < 5.0 

• -2.5 < [Ee/H] < +0.5 

However, the step sizes and ranges for log^ and [Ee/H] 
vary with effective temperature. The ASSeT spectral 
synthesis package and Kurucz model atmospheres were 


used for cool and warm stars (3,500 < Tk// < 14,000 
K), Synspec with TLUSTY model atmospheres for hot 
stars (Tpff + 15,0 00 K), and BT-Settl model spectra 
(lAllard et al. 2011|) for very cool stars (2,700 < Teff < 
+300 K). Tfiis new RV mini-grid covers a larger range of 
parameter space and with finer sampling than the first 
generation grid. In addition, a number of spectra with 
high carbon and also high a-elements are included to 
help serve as templates for carbon-rich and oxygen-rich 
stars. The synthetic spectra have a resolution of 23,500 
and are on the same logarithmically-spaced wavelength 
scale as the APOGEE combined spectra. 

We discovered that the BT-Settl spectra used for the 
coolest temperatures have a systematic, temperature- 
depende nt radial velocity o ffset on the order of 1 


km s ^ (Gottaar et al 


2014 see Eigure 7). The cause 


is not entirely clear but one possibility is a wavelength 
shift of the molecular water lines in the BT-Settl linelist. 
No RV corrections were applied in DRll+12, but future 
releases will likely have the BT-Settl spectra removed 
from the mini-grid. 


8.2.5. Synthetic Radial Velocities 

After the best fitting template is determined, each in¬ 
dividual visit spectrum is cross-correlated against this 
template to derive what we call “synthetic” radial veloc¬ 
ities. We prefer the relative velocities derived (as dis¬ 
cussed above) from the cross-correlation of each visit 
with the combined spectrum, because this should be a 
better match that does not depend on accuracy or com¬ 
pleteness of the synthetic library. However, this tech¬ 
nique can perform poorly for faint stars because the high¬ 
est S/N visit spectrum (used as the template in the first 
iteration) is still quite noisy, and especially when only a 
small number of the visit spectra are available so that the 
combined spectrum also has low S/N. In APOGEE-2, a 
significant number of faint halo stars and bulge RR Lyrae 
stars will be observed and, therefore, the RV algorithms 
will need to be improved to accomodate these faint stars, 
likely by relying more heavily on synthetic spectra. 

The synthetic RVs provide a check of the relative RVs 
for objects where there is a good library match. The scat¬ 
ter between the two types of RVs is stored in SYNTH- 
SGATTER, and when this is larger than 1 km s“^, 
the SUSPEGT_RV_GOMBINATION bit is set in the 
STARELAG bitmask. The cross-correlation functions 
from which the synthetic RVs are derived are also useful 
for detecting double-lined spectroscopic binaries (SB2s) 
because the combined spectrum will be a less helpful 
template in those cases. No automatic SB2 identifica¬ 
tion algorithm is currently in use by the pipeline, but 
the synthetic RV cross-correlation functions are saved in 
the apStar files and can be used for further inspection. 


8.2.6. B ary centric correction 

Radial velocities in APOGEE are reported with re¬ 
spect to the center of mass of the Solar System - the 
barycenter. The individual exposures are corrected for 
the relative motion of the Earth along the line-of-sight 
of the star during each observation. This is called the 
“barycentric correction” and can be calculated very accu¬ 
rately (to m s ~^ levels). Our ro utines are partially based 
on those from McCarthy (1995). When these corrections 
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Fig. 24.— Histograms of RV scatter in the three APOGEE RV 
measurements (Refined: black, Synth: red, EST: blue) for giant 
stars with multiple visits (N>3 and total S/N >20). The distribu¬ 
tions peak at '^110 m s“^ for DRDll (left) and at ~80 m s“^ for 
DR12 (right). The Refined and Synth methods have very similar 
distributions with slighter larger scatter values than EST in both 
DRll and DR12. The RV scatter is a good internal measure of 
the the internal APOGEE RV precision and long-term stability. 
Note that these values include real astrophysical variations due to 
binaries, which likely explains the long tail. 


with an rms scatter of 0.192 km s For DRll, there 
are only 15 stars in common and 

(iNidever/Chubak “ Idrii) = -0.615 ± 0.089 km 

with an rms scatter of 0.333 km s“^. Therefore, the 
changes to the RV software from DRll to DR12 im¬ 
proved both the zero-point (by ^0.25 km s“^) and dis¬ 
persion. We find no clear trends with Te// , log^f or 
[Fe/H] and do not correct the RVs for any offsets. 

8.3. RV Uncertainties 

The RV uncertainty depends on the S/N^ the reso¬ 
lution, and the information contained in the spectral 
lines themselves: A spectrum with lots of deep and thin 
lines (such as in cool and metal-rich stars) will have a 
much more accurate RV than a spectrum with few shal¬ 
low and/or wide lines (such as for metal-poor or hot 
stars). The cross-correlation RVs estimate the velocity 
uncertainty from the uncertainty in the measurement of 
the cross-correlation peak, which is partially set by the 
width of the peak. This currently systematically under¬ 
estimates the uncertainty in the RV measurements and 
will be improved in the future. 

A better estimate of the internal precision is the RV 
scatter for stars with multiple measurements. For gi¬ 
ants with a total S/N>20 and three or more visits the 
distributions peak at ^110 m s“^ for DRll and ^70 
m s“^ for DR12 (Figure |^. The EST scatters are 
slightly smaller than for tEe Refined and Synth meth¬ 
ods. This is likely because the EST method uses the 
visit spectra on their native and oversampled wavelength 
scale while both the Refined and Synth methods use the 
resampled/downsampled visit spectra. Euture improve¬ 
ments on the RV software will likely use the visit spec¬ 
tra on their native wavelength scale. The RV scatter is 
higher for (1) dwarfs, because their lines are broadened 
by rotation and higher surface gravities, and, (2) espe¬ 
cially for hotter stars with broader lines. Eignreshows 
the dependence of V^catter on nff, S/N^ and metallicity. 

Another estimate of the internal accuracy and long 
term stability of the APOGEE RVs are median plate- 
to-plate RV differences using stars in common. Eigure 
shows the histogram of median RV differences of 4317 
plate pairs with more than 50 stars in common and RV 
difference uncertainties less than 0.05 km s“^. The dis¬ 
tribution is centered around zero and has an rms scatter 
of 0.044 km s“^. This indicates that the APOGEE in¬ 
strument and the RVs have been very stable over the 
three years of survey operations. 


are applied to the absolute RVs from above we attain the 
RV with respect to the barycenter, or Vhelio foi* short. 

8.2.7. Absolute RV Zeropoint 

To ascertain the accuracy of the APOGEE velocity 
zero-point, we comp ared the APOGEE RVs t o those of 
Nidever et al. (2002) and Ghubak et al. (2012), which is 
on the Nidever velocity scale and accurate to ^30 m s“^ . 
Eor DR12, there are 41 unique stars in common between 
APOGEE and Nidever/Ghubak (7 with Nidever et al. 
and 40 with Ghubak et ah). We find 


8.4. Binarity 

The pipeline currently has no flag for binarity because 
the RV uncertainties are underestimated. In addition, gi¬ 
ants have significant astrophysical RV sc atter (“jitter' 


espec ially at the tip of the RGB (e.g., Hekker et al 


2008) that makes binary classification more complicated 


Work is ongoing to develop a reliable binary classifica¬ 
tion scheme. However, in the meantime large RV scatter 
may point to binarity. The first catalog of stellar and 
substellar companions to APOGEE stars will released in 
the near future (Troup et al. 2015, in preparation). 


(^^idever/Chubak bbR12) 0.355 dz 0.033 km S 
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Fig. 25. — The dependence of i4catter on Te.ff (left), S/N (middle) and [Fe/H] (right). Stars are selected to highlight the dependence 
in each panel. The trends are as expected with the scatter increasing for higher T^ff, lower S/N, and metal-poor stars. The blue line 
indicates median values in bins of the abscissa. 



Fig. 26.— The distribution of plate-to-plate velocity differences 
for 4317 plate pairs with 50 or more stars in common and an uncer¬ 
tainty in the velocity difference of less than 0.05 km s“^. The rms 
scatter is 0.044 km s“^ indicating that the RVs are very stable. 

All of the data are accessible through the SDSS-III 
Science Archive Server (SAS), which provides web ac¬ 
cess to the entire set of data products, ranging from 
the raw data cubes to the reduced spectra (including 
intermediate data products). These are organized in 
a directory structure and files that are described by 
the SDSS-III datamodel, which can be viewed at http: 
//data.sdssS.org/datamodel. The SAS can be ac¬ 
cessed at http://data.sdss3.org/sas, with the top 
level for APOGEE raw data at http://data. sdss3. 
org/sas/drl2/apogee/spectro/data, and the top level 
for APOGEE reduced data products at http://data. 
sdss3.org/sas/drl2/apogee/spectro/redux. Abbre¬ 
viated descriptions of the APOGEE project, instrument, 
software, and data products can be found on the SDSS- 
III DRI2 web site (http://www.sdss.org/DR12) and 
summary file descriptions can be found in the SDSS-III 
data model. 

The main products of the reduction pipeline are the 
reduced visit spectra, which are stored in apVisit files, 
and the reduced combined spectra, which are stored in 
apStar files. A webapp interface to download these for 
individual targets and bulk target lists is available at 
http: //data. sdss3. org. Otherwise, users can navigate 


their way through the SAS directory structure to find 
individual files and intermediate data products. 

Note that most of the image/spectra data products are 
multi-extension EITS files, where the contents of the ex¬ 
tensions are described in the datamodel. Eor IDL users, 
the SDSS-III apogeereduce software product contains a 
routine, APLOAD.PRO, that reads all of the extensions 
with a single command and stores the results in an IDL 
data structure. 

Parameters extracted from the spectra, along with ob¬ 
ject information, are stored in summary EITS table files. 
The allVisit file contains information about individual 
visits, while the allStar file contains information from 
the combined spectra. These tables are also loaded in 
the the Gatalog Archive Server (GAS), a database with 
a web interface at http: / / sky server. sdss3. org. 

The official APOGEE data model - defining the 
directory structure, filename convention, file formats, 
and header keywords - is available from http://data. 
sdss3.org/datamodel/files/AP0GEE_R00T/. 

10. SUMMARY 

We have described the automated APOGEE data re¬ 
duction pipeline. The basic steps are: 

1. Gollapsing the 3D up-the-ramp data cube to a 2D 
image and removing cosmic rays. 

2. Extraction of the 300 spectra from the 2D image 
and rough flux calibration. 

3. Wavelength calibration, sky and telluric correction, 
combination of dither pairs, and absolute flux cal¬ 
ibration. 

4. Combination of visit spectra and radial velocity de¬ 
termination. 

Some of the non-standard features of the pipeline are: 

• Cosmic ray rejection: Having up-the-ramp data 
gives us the ability to detect and remove most of 
the cosmic rays in our data. 

• Dither combination: The APOGEE spectra are 
slightly undersampled and, therefore, APOGEE 
observations are taken as pairs with a half pixel 
spectral dither between them. These spectra are 
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“dither-combined” with a sinc-interlace interpola¬ 
tion to create a single well-sampled spectrum. 

• Model telluric correction: We fit fairly simple tel¬ 
luric absorption models to hot star spectra across 
the plate to derive “scaling” values for each telluric 
species. Two-dimensional polynomial fits are then 
performed of the variations of these scalings across 
the plate and subsequently used to derive a model 
telluric absorption spectrum for each observed sci¬ 
ence spectrum. 

• Iterative relative RV determination: We use the 
combined observed spectrum of each star as it’s 
own RV template to derive precise relative RVs. 
This is performed in an iterative fashion since the 
relative RVs and the combined spectrum depend 
on each other. 

Euture improvements include better sky subtraction, 
persistence correction, and updates to the RV routines 
for low-S/N spectra. 
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